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Porous electrodes composed of multiphase active materials are widely used in Li-ion batteries, but their 
dynamics are poorly understood. Two-phase models are largely empirical, and no models exist for three 
or more phases. Using a modified porous electrode theory based on non-equilibrium thermodynamics, 
we show that experimental phase behavior can be accurately predicted from free energy models, without 
artificially placing phase boundaries or fitting the open circuit voltage. First, we simulate lithium interca¬ 
lation in porous iron phosphate, a popular two-phase cathode, and show that the zero-current voltage gap, 
sloping voltage plateau and under-estimated exchange currents all result from size-dependent nucleation 
and mosaic instability. Next, we simulate porous graphite, the standard anode with three stable phases, 
and reproduce experimentally observed fronts of color-changing phase transformations. These results 
provide a framework for physics-based design and control for electrochemical systems with complex 
thermodynamics. 

© 2014 Elsevier Ltd. All rights reserved. 


1. Introduction 

In order to develop safer, longer-lasting, and more energy-dense 
batteries, it is crucial to understand their thermodynamic behavior 
out of equilibrium. Many battery materials exhibit multiple phases 
with varying composition, voltage, and temperature [1-4], driven 
by electrochemical reactions [5], In Li-ion batteries, complex phase 
transformations are triggered by lithium intercalation reactions. 
The standard anode material, graphite, passes through three or 
more phases [6] with observable color changes [7], while popular 
cathode materials, such as olivine phosphates [8-11] and transition 
metal oxides [12,13], exhibit two-phase separation across single 
particles [14,15] and porous electrodes [16], Conversion reactions 
can also lead to complex phase behavior in lithium-sulfur batteries 
[17-19] and lithium-oxygen batteries [20,21], 

Phase transformations pose a major challenge for mathematical 
models used to design, characterize, and control Li-ion batteries 
[22], Systems-level models neglect phase transformations alto¬ 
gether and rely on empirical constructs, such as current-dependent 
particle sizes [23,24], Classical porous electrode theory (PET) cap¬ 
tures the microscopic physics of diffusion and reactions [25-27], 
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but does not consistently describe multiphase thermodynamics 
[28,29,5,30], In all cases, the voltage plateau for two coexisting 
phases is fitted to a unique (single-phase) function of the state of 
charge. For two-phase materials, such as Li x FeP0 4 (LFP,X~0,1) [8], 
spherical “shrinking core” [25,31 ] or planar [32-34] phase bound¬ 
aries are imposed in the active particles, but no such models are 
available for materials, such as graphite, with three or more phases. 

Recent work on LFP has shown that the voltage plateau is 
an emergent property of the active particles [11,35-41] and the 
porous electrode [42,43,30,44,45] for any material with multi¬ 
ple stable phases at different compositions. In single particles, 
phase separation occurs within the miscibility gap (range of 
unstable homogeneous compositions), which depends on temper¬ 
ature and particle size [11,37], overpotential [46,47,11,48], current 
density [36,37,49,50], and particle connectivity [51 ]. In porous elec¬ 
trodes, a collection of particles entering the spinodal gap (linearly 
unstable compositions) can also undergo discrete transformations 
[42,43,16,52], but the roles of nucleation [48,44], ion transport 
[30,53], and heterogeneity [44,16] are just beginning to be under¬ 
stood. 

In this work, we present a predictive theory of phase transforma¬ 
tions in Li-ion batteries. The mathematical framework, combining 
classical porous electrode theory [27] with non-equilibrium ther¬ 
modynamics [5], is described in a companion paper [30], This 
new approach effectively homogenizes microstructural simula¬ 
tions [54,55], even with phase separation [56,51 ], at a tiny fraction 
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of the computational cost. Here, we introduce simple free energy 
models for LFP (two phases with coherent nucleation) and graphite 
(three phases, neglecting nucleation) and use our modified porous 
electrode theory (MPET) to predict representative experimental 
data [42,7], 

2. Iron Phosphate: Two Phases 

The understanding of binary phase separation in single particles 
is rapidly advancing [5 ]. Phase-field models of LFP have been 

developed for isotropic spherical particles [47,11,57,40,41], which 
may be relevant for large particles or agglomerates, while new 
modes of intercalation have been identified [5,39] for anisotropic 
single-crystal nanoparticles [10,58], At low currents, intercalation 
waves sweep across the active (010) surface [59,60,36,61 ] and relax 
to striped patterns along {100} [14] or {101} [15] planes due to 
coherency strain [37], Once the local current density exceeds the 
exchange rate, phase separation is suppressed [36,37], A “single¬ 
phase transition path” has also been suggested based on the 
large barrier for bulk nucleation [62], but there is new evidence 
(strengthened below) that two-phase nucleation does occur, via 
the instability of surface layers [48], The nucleation barrier is size- 
dependent and vanishes below a critical particle size around 22nm 
[48], consistent with experiments [63], 

Much less is known about binary porous electrodes. The semi¬ 
nal work of Dreyer et al. demonstrated a zero-current gap between 
voltage plateaus for charging and discharging [42], which they 
attributed to discrete spinodal decompositions among bi-stable 
homogeneous particles, analogous to filing a balloon array [43], 
Independently, Burch [52] observed the “mosaic instability” in 
simulations of a collection of phase-separating particles described 
by the Cahn-Hilliard-reaction model [64,5,40,41 ], which exchange 
ions through an electrolyte reservoir at constant total current. (See 
Ch. 9 of Ref. [52].) Recently, Bai and Tian (BT) interpreted cur¬ 
rent transients in terms of population dynamics for nucleation 
and growth in a set of discrete particles [44] in place of the 
Kolmogorov-Johnnson-Mehl-Avrami statistical theory for a con¬ 
tinuous electrode [65], The BT theory has already led to accurate 
measurements of mechanical deformation [45] and charge-transfer 
kinetics [66], and the mosaic picture is directly supported by the 
first porous-electrode imaging experiments of Chueh et al. [16]. 
The data reveal mostly single-phase particles (x**0, 1) and some 
two-phase particles with planar (not core-shell) phase boundaries, 
consistent with theory [59,36,37,48], Macroscopic phase gradients 
are also observed [16], but have never been modeled. 

Here, we use MPET [30,5,53] to re-interpret the data of Dreyer 
et al. [42] at very low C-rates, C/n, discharging the capacity in n = 
50-1000 hours. The porous cathode is partitioned into finite vol¬ 
umes, each containing a representative LFP particle of random size 
(log-normally distributed) and a realistic shape (C3) [14,67,48], 
The separator and lithium anode (at constant potential) are also 
modeled. The homogeneous free energy 

g(x) = k B T [xlnx + (1 -x)ln(l -x)]±£2x(l -x) (1) 

and diffusional chemical potential 

^ = S = kBTln ( T3^)+^ (1 - 2x) ( 2 ) 

(per site) describe a regular solution of particles and vacancies with 
mean interaction energy £2, previously fitted to Li-solubility data 
versus temperature and particle size [37], 

Based on porous electrode imaging [16] and short diffusion 
times (~ ms) in nanoparticles [68], we assume fast single-particle 
transformations compared to the C-rate and do not resolve the 
internal dynamics of each particle. Instead, we replace each par¬ 
ticle with an effective homogeneous solid solution (the “pseudo 


capacitor approximation” [30]), whose electrochemical response 
is tuned to the results of a realistic phase-field model [48], In par¬ 
ticular, the regular solution parameter £2 is varied with particle size 
to match the spinodal to the coherent nucleation voltage for dis¬ 
charging^ insertion) [48], predicted from ab initio surface energies 
[69] and elastic constants [70], (See Appendix A.) As a first approx¬ 
imation, the same size-dependent nucleation voltage is also used 
for charging. The local Nernst equilibrium voltage (relative to the 
lithium metal anode) is defined by the chemical potential, Ve Q = 
V® — fi/e, relative to a formal reference voltage V® at x = 1 /2 in the 
middle of the plateau [5,30,36], The local overpotential r) = V- V eq 
determines the reaction rate via generalized Butler-Volmer kinetics 
from the local interfacial voltage Atf> [5], Ion transport is governed 
by a standard Nernst-Planck electrolyte model [30,27], The sim¬ 
ulations follow the experimental charge/discharge protocol, but 
a wider range of electrode filling fraction is covered in order to 
minimize startup effects and reach the steady voltage plateau well 
before comparison with the experimental data (Fig. la). 

Our MPET has only three fitting parameters: the mean 
(fi = 28nm) and standard deviation (a = 3.5nm) in a log-normal dis¬ 
tribution of the cycled particle size (defined by the short axis of the 
C3 shape in the [100] direction, which is roughly half of the long 
axis length) and a series resistance (R s = 3.9£2g). Unlike traditional 
models, the voltage plateau is not fitted, but predicted, and phase 
transformations occur spontaneously. Material properties of LFP 
were previously determined by ab initio calculations and exper¬ 
iments [37,48], except for the exchange current (see below). All 
porous electrode properties are known or estimated (electrolyte 
diffusivity, active material loading, porosity, thickness, tortuosity), 
but do not significantly affect the results at low rates. 

With so few adjustable parameters, the agreement between the 
theory and experiment in Fig. 1(a) is remarkable. The discharge 
plateaus are reproduced with millivolt accuracy, including a slight 
tilt that had escaped notice. The charging plateaus are also well 
described, considering that a separate model was not developed 
for charging nucleation. (The charging data also strangely overlaps 
for C/200 and C/131, separate from C/1000.) 

An unexpected finding is the effect of particle-size variance on 
the battery voltage. For a single particle size, the voltage remains 
constant during phase transformation at zero current. In a realistic 
heterogeneous composite, the voltage plateau tilts as particles fill in 
order of increasing nucleation overpotential [48], from smallest to 
largest. Consequently, the zero-current voltage profile can be used 
to infer the particle size distribution, either by simulations ( Fig. 1 ) or 
from a simple analytical formula derived in Appendix C. As shown 
in Fig. 1(b), the theory predicts that only the smallest of the active 
particles were cycled in these experiments, having (i = 28 nm and 
a = 3.5 nm. The reported range of 50 to 1000 nm [42] may have 
reflected agglomerates or limited imaging resolution, since the 
smaller values also lead to accurate predictions of the voltage gap 
between charge and discharge versus current, as shown in Fig. 1(c). 
Additional support for our interpretation comes from recent analy¬ 
sis of the nucleation voltage versus LFP nanoparticle size (Fig. 4c of 
Ref. [48]), where extensive literature data mostly collapses with¬ 
out any parameter fitting. The present data for a partially cycled 
electrode is anomalous using the reported mean particle size, 
but falls back on the universal scaling curve using the size range 
28 ±3.5 nm inferred here from the tilt of the voltage plateau. The 
theory of size-dependent coherent nucleation also explains why 
the observed zero-current gap of 20 mV is much smaller than 
the bulk spinodal gap of 74 mV inferred from thermodynamic 
data [37], 

Another unexpected finding is the effect of particle-size vari¬ 
ance on the spatial profile of the phase transformation, shown in 
Fig. 1(d). With identical particles (cr = 0), the mosaic instability is 
localized in a thin reaction front propagating from the separator, 
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Fig. 1. Simulations with modified porous electrode theory for two-phase Li„FeP 04 , compared to the experimental data of Dreyer et al. [42], (a) Charge-discharge cycles at 
C/1000, C/200, and C/131, (b) Predicted effect of the particle-size standard deviation a (in a log-normal distribution) on the discharge curve at C/200 for the same mean 
size ft = 28nm. The simulation results are smoothed; the raw data appears in Fig. 4. (c) Charge-discharge voltage gap, with additional C/10 and C/50 experimental data. 
Theoretical curves are shown for a single particle size (cr=0) and for the fitted size distribution in the inset, (d) Simulated lithium profiles across the porous cathode at 50% 
filling for the cases in (b), showing the transition from a sharp reaction front at a - 0 to homogeneous filling in order of increasing particle size with increasing a. A video of 
the full C/200 discharge simulation from (d) synchronized with the voltage profile in (b) for <r=3.5 nm is linked from the online version of the article and also available at 
https://www.youtube.com/watch7v~PlE31ZaJXQ. 


whose width is proportional to the current, and voltage fluctua¬ 
tions are enhanced by blocks of transforming particles [30,5], With 
non-identical particles, a small standard deviation (as small as a = 1 
nm, not shown) suffices to spread the phase transformation over 
the entire electrode and smooth the voltage profile. The particles 
fill in spatially extended groups according to size, starting with the 
smallest, with a gentle gradient in concentration away from the 
separator. This theoretical picture is consistent with the imaging 
experiments of Chueh et al. [16], 

These findings can drastically alter the interpretation of electro¬ 
chemical data for multiphase porous electrodes. Classical methods 
[71 ], such as intermittent titration or impedance spectroscopy [72], 
are based on reaction-diffusion models that predict spatially uni¬ 
form reactions at low current. For solid-solution materials, such 
as silicon nanowires, this assumption is justified and leads to very 
accurate impedance analysis [73], but for phase-separating mate¬ 
rials, mosaic instability implies that low currents may only probe 
small fractions of the total active area and particle-size distribution. 
The reaction rate can be grossly under-estimated if the full internal 
surface area is assumed to be active. Indeed, classical PET yields very 
small exchange current densities (3xl0 -6 [25] and 5.4xl0 -5 [31] 
A/m 2 ), while our MPET uses 7x10 -3 A/m 2 , which seems more rea¬ 
sonable for a high-rate cathode. Much larger values >10 A/m 2 have 
also been reported [32,33] that are >10 7 times the PET values, so it is 
clearly important to account for the mosaic instability using MPET. 
Since the particle-size variance promotes uniform electrode inter¬ 
calation, our results show that the simple BT population-dynamics 


model [44] can also be a good approximation at low rates and/or 
early times (without electrolyte depletion) to replace the Cottrell 
equation [ 71 ] in chronoamperometry as a method of measuring the 
exchange current, e.g. 10 -4 A/m 2 in Ref. [66], 

It is important to stress that our successful fitting of MPET is 
for low-rate experiments, probing the slow dynamics of electrode 
phase transformations. At high rates, the original MPET formula¬ 
tion with Butler-Volmer kinetics and no statistical hetereogeneities 
[30] can have difficulty fitting experimental data [53], The frame¬ 
work of MPET is very general, however, and it is likely that some 
more effects must be added to describe the full range of operat¬ 
ing conditions. This work shows that a distribution of particle sizes 
should be considered, while other studies show the importance of 
electron transport limitation [74] and higher-resistance Marcus- 
Hush-Chidsey kinetics of electron transfer from the carbon coating 
of LFP nanoparticles [66], 

3. Graphite: Three Phases 

Traditional battery models cannot describe active materials 
with three or more phases, and yet many insertion compounds have 
this property. Empirical approaches, such as “shrinking annuli” [75] 
that generalize the two-phase “shrinking core", are difficult to jus¬ 
tify theoretically and require inconsistently fitting the open circuit 
voltage “staircase" (set of plateaus for multiple phase transforma¬ 
tions) as a unique function of the mean concentration. On the other 
hand, our MPET has no such limitations and simply requires fitting a 
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Fig. 2. Graphite MPET simulations vs. experiment, (a) Open circuit voltage “staircase” for lithium insertion in graphite (vs. Li metal) [6], reproduced by the MPET model (solid 
curve), whose homogeneous free energy is the dashed curve, (b) Contour plot of the free energy model versus the filling fractions of two adjacent, periodically repeating 
layers with repulsive interactions across layers. The model has local minima near (0,0) for the empty phase (blue), (1,1) for stage I (gold) and (1,0) or (0,1) for stage II 
(red). These minima are connected by common tangent constructions for two-phase coexistence, which force the system from empty to stage II to stage I during filling, 
(c) Simulated red/gold interface versus time (solid curve) in the experiments of Harris et al. [7], compared to the experimental data, (d) Experimental image (above) and 
simulated color profile (below) at the same moment in time, capturing the blue/red and red/gold interface positions and noise. A video of the full simulation synchronized 
with the experiment is linked from the online version of the article and also available at https://www.youtube.com/watch?v=vroBPmOfODk. 


homogeneous free energy model with three or more local minima 
to the phase diagram; the voltage staircase is then an emergent, 
rate-dependent property of the material. 

To illustrate this new approach, we apply MPET to the impor¬ 
tant case of graphite, a complex intercalation compound [76] that 
has at least five distinct phases during lithium insertion, Li x Cg, as x 
varies from 0 to 1 [6], The open circuit voltage (Fig. 2(a)) exhibits 
two wide plateaus from 1/4 tox«* 1/2 (stage II) tox« 1 (stage I), 
corresponding to the stable states with every third, second, and sin¬ 
gle layer between graphene planes filled with lithium, respectively. 
Higher stages with 0 < x < 1 /4 involve correlated interlayer domains 
[77] that are harder to discern experimentally and less important 
for intercalation dynamics. As a first approximation, therefore, we 
consider only the three phases at x«*0, x«*l/2 and x«*l, using a 
simple free-energy model proposed by Bazant [78], 

In order to capture two voltage plateaus between these states, 
the model introduces two order parameters, (xi, x 2 ), that represent 
the lithium fractions in a pair of periodically repeated layers, where 
x = (xi +x 2 )/2. The inset of Fig. 2(a) illustrates this physical picture 
of the three phases at the endpoints of the voltage plateaus. The 
free energy per site in the layer pair, 

g(x t ,x 2 ) = «(xi) + g(x 2 ) + g int (xi ,x 2 ), (3) 

has local minima near (1,1) for stage I (x r* 1) and (0,0) for the low 
density phase (x?»0), as well as equivalent minima near (1,0) and 
(0,1) for stage II (x«* 1/2), as shown in Fig. 2(b). Each layer i = 1, 2 
has a double-welled homogeneous free energy, g(Xj), with minima 


near x,- = 0,1 favoring full or empty states. Stage II with every other 
layer full is stabilized by a repulsive interaction, g, n t(xi, x 2 ), between 
adjacent layers. 

Building on the LFP model above, each layer is treated as a reg¬ 
ular solution, 

I(Xf) = k B T [X; lnx; + (1 - x f )ln(l - x,)] + S2 a x,(l -x,) (4) 

where £2 a > 0 controls the width of the each voltage plateau. The 
interaction energy has a similar polynomial form 

& nt (xi,x 2 ) = S2 b x 1 x 2 + S2 c xi(l -xi)x 2 (l -x 2 ) (5) 

where £2j, > 0 is a repulsion energy between two particles at the 
same site in adjacent layers, which sets the different between the 
two voltage plateaus. The second term with S2 C >0 (not in the 
original model [78]) is an additional repulsion between adjacent 
particle-vacancy dipoles, which penalizes partially filled layers and 
further stabilizes stage II with Xi =/= x 2 . Intercalation reactions are 
allowed to proceed into each layer independently, as if it were 
a separate reactant. Each layer thus has its own local overpoten¬ 
tial, iji = v-v eq ,u and Nernst voltage, V eq j = V® - /t;/e, defined by its 
diffusional chemical potential, 

Mi = || = /If + O* + fi c x,(l - XjX 1 - 2Xj) (t#l) (6) 

where 

Hi = || = k B nn (y^t) + n«(i - 2 x,). 


(7) 



















T.R. Ferguson, M.Z. Bazant/ Electrochimica Acta 146 (2014) 89-97 


93 


The free-energy model has four parameters, £2 a , £2 C and 

V®, that are fitted to reproduce the open circuit voltage at a given 
temperature (Fig. 2(a)). 

Graphite is an electrochromic material that changes color upon 
lithiation [79-81 ], making it possible to visualize its phase transfor¬ 
mations. The low density state is black and switches to blue from 
a small value of x<0.05 up to 1/4. Stage II is red and extends 
just past x?»l/2, while stage I is gold and covers the widest volt¬ 
age plateau up to x = 1. These colors are indicated by circles around 
the corresponding minima of the model free energy in Fig. 2(b). 
Slow lithium insertion follows the path of lowest free energy of 
the convex hull, including common tangent planes between the 
local minima that represent two-phase co-existence. Lithiation 
starts near the blue circle (homogeneous empty state) and pro¬ 
gresses towards one of the minima inside the red circles (replacing 
the empty state with stage II) then towards the minimum inside 
the gold circle (replacing stage II with stage I). The contour plot 
also shows the tilt of the free energy along the (1,1) direction, 
which is controlled by and leads to the different voltage plateau 
values. 

The experiments of Harris et al. [7] visualizing coloration dur¬ 
ing lithium intercalation in an “unrolled” porous graphite electrode 
provide a unique opportunity to test our MPET model for a material 
with more than two phases, for the first time. Beautiful experimen¬ 
tal movies (posted at http://lithiumbatteryresearch.com) show the 
nucleation of a sequence of thin reaction fronts propagating diffu¬ 
sively from the separator into the porous electrode during lithium 
insertion, switching the color of the graphite particles from black to 
blue to red to gold in discrete, stochastic transformations (Fig. 2)(d). 
The battery is discharged at a constant potential of 2 mV, very close 
to short circuit, which leads to a large initial flux of lithium into the 
particles. 

Since the characteristic time for diffusion across this long elec¬ 
trode (~1 mm) is on the order of hours, the electrolyte quickly 
becomes depleted of salt. The propagation of the reaction front 
is thus limited by the lithium diffusion that feeds it (from both 
sides), leading to the same square root of time scaling of the front 
position and diffusion layer width [7], as in diffusion-limited corro¬ 
sion of porous electrodes [82,83]. Relative to electrolyte diffusion, 
reactions and solid diffusion are very fast, so the pseudo-capacitor 
approximation can be safely made in our MPET [30], as in the case 
of LFP nanoparticles above. Ohmic losses for electrons can also be 
neglected since the porous graphite layer is very thin and sits on 
the current collector. 

The MPET simulations also include the separator and lithium 
metal anode and apply the same constant potential (2mV) as in the 
experiments. In order to compare with the experimental images, 
the solid lithium concentration is converted to three colors, con¬ 
sistent with the free energy model and experimental observations: 
blue 0 < x < 0.3, red 0.3 < x < 0.6, and gold 0.6 < x < 1. See Appendix 
B for further details. 

With essentially no adjustable parameters, the agreement 
between theory and experiment in Fig. 2(c)-(d) is remarkable. Once 
the free energy model is determined by the open circuit voltage, 
the only fitting parameter is the effective diffusivity D e jj of the 
electrolyte (since both reactions and solid diffusion are fast), but 
its value must remain close to theoretical estimates. Assuming 
the Bruggeman relation, D e g = ^J 2 D, with an ambipolar diffusiv¬ 
ity D = 4.6x 10 -10 m 2 /s (consistent with experimental values for 
ethylene carbonate / dimethyl carbonate solution [84]), an excel¬ 
lent fit of the position of the red/gold interface versus time (Fig. 2(c)) 
is obtained with a reasonable porosity, € p = 0.4. As anticipated by 
the experimentalists [7], the porosity is much larger than in com¬ 
mercial graphite porous electrodes (15%) because the electrode 
consists of roughly a monolayer of irregular 5-20/xm particles [7] 
without any additives compressed between parallel flat surfaces 


for visualization purposes. The porosity of random sphere pack¬ 
ings (measured by Voronoi tessellation) is known to increase from 
37% in the bulk to around 50% when pressed against a flat surface 
[85], and since irregular, polydisperse particles pack more densely 
than spheres [86], the fitted value of 40% is quite consistent with 
the experimental geometry. 

Fitting the diffusive motion of a single interface may not seem 
so surprising (and can also be done by an ad hoc diffusion equation 
for coloration [7]), but our MPET accurately predicts the nucleation 
and propagation of two stochastic reaction fronts, the red/gold and 
blue/red interfaces, over the entire recorded time without fitting 
any additional parameters. (See the supporting movie and Fig. 2(d).) 
Moreover, a close examination of the reaction fronts shows that 
the mosaic instability, i.e. stochastic filling of discrete particles 
[42,52,44,30], is quite similar in both movies, from experiment and 
simulation. 

4. Conclusion 

In this paper, we show that MPET [30] based on electrochemical 
non-equilibrium thermodynamics [5] is able to accurately simulate 
two fundamental experiments with multiphase porous electrodes 
[43,7 ] that traditional porous electrode theories could not describe. 
The advantage of MPET is that it couples the thermodynamics of the 
active material to electrochemical transport and reaction kinetics. 
Complex dynamical phenomena, such as nucleation, phase growth, 
mosaic instability, and voltage hysteresis, are then predicted by the 
model, rather than artificially imposed on the system. 

The fundamental input to MPET is a free energy model for the 
active material, inferred from the phase diagram and open cir¬ 
cuit voltage. The cell voltage is predicted as an emergent property 
of the porous electrode, rather than fitted to experimental data 
as an effective property of the active material, as if it remained 
homogeneous and never phase separated. By properly accounting 
for non-equilibrium thermodynamics, MPET provides a promising 
framework for battery modeling to optimize the cell design, predict 
and control performance, monitor the internal state, and improve 
safety under diverse operating conditions. 

Video Content 

Movies of the simulations in Figures 1 and 2 are linked 
from the online version of the article and can also be found 
at the Battery Modeling Group YouTube page, located at 
https://www.youtube.com/user/BatteryModelingGroup. 
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Appendix A: Iron Phosphate Simulations 

The general simulation framework of MPET is described in 
a companion paper [30] and is beginning to be used by other 
researchers [53], In this Appendix, we give some of the parameters 
and simulation details not covered in the main text. 

Cogswell demonstrated that single particle voltage profiles are 
always tilted due to coherency strain and that the overshoot from 
the standard potential on discharge varies with particle size. This 
overshoot, which can be thought of as the half gap (i.e. the overpo¬ 
tential required to drive lithiation), depends on particle size since 


94 


T.R. Ferguson, M.Z. Bazant / Electrochimica Acta 746 (2014) 89-97 



Wetted Area/Volume (nm 1 ) 



Fig. 3. Size-dependent regular solution model to approximate coherent nudeation 
in nanoparticles, (a) Theoretical prediction of the critical nudeation voltage in LFP 
(blue), decreasing linearly with area to volume ratio [48], and the effective size- 
dependent regular solution parameter Q eff . (b) Equilibrium voltage profiles for 
different particle sizes from the size-dependent regular solution model fitted in 
(a). 


the surface wetting represents a larger percentage of the total vol¬ 
ume for smaller particles. 

As noted in the main text, the particles are modeled as effec¬ 
tively homogeneous (“pseudo capacitor approximation” [30]), but 
with an equilibrium voltage profile that approximates coherency 
phase separation within the particle [37], We simply adjust the reg¬ 
ular solution parameter with particle size to place the unstable 
spinodal points (extrema of the chemical potential) at the size- 
dependent nudeation voltage predicted by Cogswell and Bazant 
[48] (Eq. 9 below), as shown in Fig. 3. The value of £2 can be 
determined from an exact expression for the voltage gap between 
spinodal points in the regular solution model, 


AVgap = 


2k B T 

e 




Q? - 20. - 2tanh -1 



(8) 


In this approach, we neglect asymmetry between charge and 
discharge in the single-particle properties. 

For purposes of calculating size-dependent nudeation barriers 
[48], the size of each representative particle in a finite volume of 
the porous electrode is sampled from a log-normal distribution. 
The particle shapes are all assumed to be the C3 shape [67,48], The 
particles are plate-like and all assumed to be 20 nm thick in the 
[010] direction [67,87], The wetted surface area to volume ratio is 
calculated via A/\7=3.6338/L, where L is the size of the particle in 
the [100] direction. [48], 



ent particle size distributions without any smoothing of the data. In order to convey 
the effect of particle size and extract the voltage gaps between charge and discharge, 
the same data is smoothed in Fig. 2(b). 


Given the very slow charge/discharge rates in the experiments 
[42], the electrolyte properties are not very important since the 
electrolyte does not deplete. For completeness, however, suitable 
numbers are chosen for transport properties. An ambipolar diffu- 
sivity of 1.5x10 10 m 2 /s is used for the 1M LiPF 6 in EC/DEC electrolyte. 
[84], A transference number does not seem to be available for this 
electrolyte, so a value of 0.35 is assumed, consistent with other 
typical battery electrolytes. [1] The electrode is assumed to be 50 
/cm long with a 25 /cm separator. The volume fraction of the active 
material is assumed to be 0.5, and the porosity is 0.4. 

The simulation is run by starting with a fully charged electrode, 
which is then discharged to a filling fraction of 0.2 at a C/10 rate. 
The electrode is then relaxed by simulating a zero current for a long 
period of time (roughly 3.5 hours). The electrode is then discharged 
to a filling fraction of 0.7, relaxed, then charged back to 0.2. When 
calculating the voltage gaps, the simulation results are smoothed 
to make the value consistent, to account for averaging over a large 
number of particles of different sizes that cannot be captured in our 
simulations with a relatively small set of discrete particle sizes. The 
simulation results from Fig. 2(b) without smoothing are shown in 
Fig. 4. A constant contact resistance of 3.9 £2g is inferred by fitting 
the data, as one of only three fitting parameters, along with the 
mean and variance of the particle size in a log-normal distribution. 

Appendix B: Graphite Simulations 

The graphite free energy model of Bazant [78] described in 
the main text is adjusted so that a slow C/1000 MPET simulation 
(without transport limitation) fits open circuit voltage staircase 
(Fig. 2(a)), yielding the parameters £2 a = 3.4kBT, Qf, = 1.4kgT, and 
V° = 0.1366 V. The fourth parameter, S2 c = 301<bT, assigns an extra 
energy to the homogeneous state that serves to push the parti¬ 
cle toward stage II at intermediate filling fractions, but its precise 
value is relatively unimportant here. Each representative layer is 
modeled as having its own reaction rate, determined by generalized 
Butler-Volmer kinetics [5], The slow discharge simulation (which 
allows the transport effects to be neglected) is used to fit the model 
parameters. 

Once the free energy parameters are obtained, a constant poten¬ 
tial discharge simulation at V= 2mV is run for a full cell simulation 
with a lithium metal anode and a graphite cathode to simulate the 
experiments of Harris et al. [7], Transport effects on the lithium 
metal electrode are assumed to be negligible, and the porosity on 
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that side is treated as unity to model the free electrolyte chan¬ 
nel above the unrolled current collector. The important transport 
effects are on the graphite side, where there is sharp electrolyte 
depletion from the separator to the intercalation front. It is nec¬ 
essary to simulate a full electrode in order to obtain a realistic 
diffusivity from fitting. In a half cell neglecting the anode, large 
salt concentration gradients form across the separator and lead to 
artificially small fitted diffusivities. 

The separator thickness from images provided in the original 
paper is estimated to be 1.23 mm. [7] The electrodes are assumed 
to be on the order of 1.2 cm (or 10 times the separator thickness). 
The total length is not important since the coloration dynamics are 
observed in the first couple millimeters. The diffusive time across 
the separator is on the order of one hour, which is rate limiting com¬ 
pared to the much faster reactions. The exchange current densities 
for both graphite and lithium metal are arbitarily set the inverse 
diffusion time across the separator, which corresponds to approxi¬ 
mately 1.4 A/m 2 for 5 /cm spherical particles. (Particle sizes of 5-20 
/cm are observed). The volume fraction of graphite in the electrode 
is assumed to be 0.8, and the Bruggeman relation is used to model 
porous transport effects [30], 


Appendix C: Theory of the Tilted Voltage Plateau 

Our MPET model reveals an unexpected effect of statistical vari¬ 
ations in the active particle size on the open circuit voltage of a 
battery, namely a slight tilt of the voltage plateau between two sta¬ 
ble phases. This feature is present in published data for LFP at very 
low rates [42], but has apparently gone undetected until this work 
(Fig. 1 ). In this section, we provide a simple analytical theory of the 
voltage profile versus mean filling, V(x), near open circuit condi¬ 
tions, based on the particle size distribution and the recent theory 
of size-dependent nucleation [48], 

Let A0 = V® - V- IR S be the voltage drop relative to the formal 
reference potential of two-phase coexistence V®, accounting for 
a small Ohmic loss from the overall series resistance R s . For low 
rates, close to open circuit conditions, we neglect any concentration 
polarization within the electrode. Let n(L) be the probability density 
function for the linear size L of an active particle. 


Following Cogswell and Bazant [48], we assume that nucleation 
occurs on certain surfaces wetted by the new phase, e.g. the high 
density phase of LFP during battery discharge. The critical voltage 
for nucleation, A0*, which corresponds to the coherent miscibility 
limit, decreases linearly with the wetted area to volume ratio, L//8, 

A</>* = At/,% (l - |0 (9) 

where fi is a constant for a given particle shape (the volume to 
wetted area ratio per linear size, assumed to be the same for all 
particles), A0*oo is the critical voltage in the large-size limit (set 
by elastic strain energy), and L c is the critical particle size below 
which the nucleation barrier vanishes (due to the dominance of 
surface energy). Simple analytical formulae are available to predict 
A iplc and L c from the fundamental thermodynamic properties of 
the active material (elastic constants, misfit strain, surface ener¬ 
gies, interfacial tension, miscibility limit, etc.), and in the case of 
LFP, the predicted values A0J, = 37 mV and L c = 22 nm lead to a 
remarkable data collapse of reported nucleation voltages for differ¬ 
ent mean particle sizes. For the C3 shape of LFP particles [67,48 ], the 
geometrical parameter is f)= (3.6338) -1 = 0.2752, if L is the particle 
length in the [100] direction. 

At a given voltage drop in the two-phase region, 0 < V < A0^, 
near open circuit conditions, all particles with I<L*(A0) are filled, 
and the others empty, where 


I*(A0) = 


pic 

1 - A0/A0*, 


(10) 


is the size of the particles at the nucleation voltage, L = L*, undergo¬ 
ing filling transformations. The mean filling fraction of the electrode 
is then given by 


f^niLpdL 
Io° n(L)L 3 dL 


( 11 ) 


which is an implicit formula for the “tilted voltage plateau” versus 
state of charge, V(x). The voltage profile is generally nonlinear, but 
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its slope at a given point (inverse of the pseudo-capacitance) is 
given by 

dA0 v P (l- A0/A<fe) 5 Afe, (12) 

* n(L*(A0))(j8L c ) 4 

where the third moment of the particle size distribution, v p = 
J 0 °° n(L)L 3 dI, is an effective mean particle volume. 

Figure 5 shows the effect of the shape of the particle size dis¬ 
tribution on the profile of the open circuit voltage versus state of 
charge. A unimodal distribution, such as a Gaussian or log-normal, 
leads to a tilted voltage plateau with a slightly nonlinear profile, as 
shown in Fig. 2, but a bimodal distribution leads to a voltage stair¬ 
case with two tilted plateaus close to the mean voltages where the 
two different particle sizes transform (at the size-dependent coher¬ 
ent miscibility limit). Interestingly, this statistical mechanism for 
a voltage staircase is distinct from that of the graphite discussed 
in the main text, which has to do with transitions between three 
or more stable phases. In principle, a very precise measurement 
of the tilted voltage plateau can be used to infer the particle size 
distribution, by solving an inverse problem given by Eq. (11). 
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